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Abstract. We apply the Bethe-Peierls approximation to the problem of the inverse 
Ising model and show how the linear response relation leads to a simple method to 
reconstruct couplings and fields of the Ising model. This reconstruction is exact on 
tree graphs, yet its computational expense is comparable to other mean-field methods. 
We compare the performance of this method to the independent-pair, naive mean- 
field, Thouless- Anderson-Palmer approximations, the Sessak-Monasson expansion, and 
susceptibility propagation in the Cayley tree, SK-model and random graph with 
fixed connectivity. At low temperatures, Bethe reconstruction outperforms all these 
methods, while at high temperatures it is comparable to the best method available 
so far (Sessak-Monasson). The relationship between Bethe reconstruction and other 
mean-field methods is discussed. 
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The enormous and ongoing growth of data in molecular biology has led to a surge 
of interest in inverse problems. Examples are the reconstruction of gene regulatory 
networks from expression levels [1], the identification of neural interactions from 
neural spike recordings [2], and the prediction of protein-protein interactions from the 
evolutionary correlation between amino acids [3]. The paradigm of such inverse problems 
is the inverse Ising model, where the parameters of the Ising model Hamiltonian 
(couplings and fields) are to be inferred from a set of spin configurations drawn 
independently from the equilibrium measure. The goal is to use correlations observed 
between spins to infer the underlying, unknown couplings and fields. This problem is 
intrinsically hard, as pairs of spins can exhibit correlations without interacting directly 
with each other. 

In this note, we consider the inverse Ising model at the level of the Bethe-Peierls 
approximation (BP) and show how the linear response approach [4, 5, 6] leads to a 
reconstruction of the Ising model that is efficient, straightforward and outperforms 
currently available mean-field-like methods in benchmarks for strong couplings (and 
does as well as they do at weak couplings). 

Consider the Boltzmann measure 

p[s] = ie-^W (1) 

of the Ising model with H[s] = hiSi — JijSiSj, specifying the statistical weight 
of a configuration s = (si, S2, sat) given couplings Jij and local fields hi (at inverse 
temperature (3 = 1), where Z is the partition function. The indices (ij) run over all 
pairwise interactions. The inverse Ising problem is then to recover the couplings and 
fields from a given set of M spin configurations D = {s^, s^, s*^} drawn independently 
from (1). The log- likelihood of couplings and fields given such a set of observations, 
normalized by M, is 

L{{sn\{h., j^A) = - In ^[{/^., j.,}]+E h^j^ E ^r+E J^^l^ E ^r^i -(2) 

i M (ij) A* 

Maximizing the log-likelihood with respect to the Ising model parameters Jij and hi 
leads to 

niiiihi, Jij}) = mf , 

Cij{{hi, Jij}) = C^j , (3) 

where = (sj) and Cij = (siSj) — {si){sj) are the magnetizations and connected 
correlations under the Boltzmann distribution (1). The right-hand sides are the 
corresponding quantities estimated from data, mf = J2fj, -sf and Cj^ = J2fj, sfsj — 
mfmf [7]. 

Many different approaches to find the couplings Jij and fields hi by maximizing 
the log-likelihood (2) or by directly solving the self- consistent equations (3) have 
been taken, including gradient descent with Monte-Carlo simulation [7], independent- 
pair approximation (IP) [8], naive mean-field (MF) [4], Thouless- Anderson-Palmer 
approximation (TAP) [4, 9], Sessak-Monasson expansion (SM) [10], susceptibility 
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propagation (SusP) [11, 12] and more recently an adaptive cluster expansion [13] 
and pseudo-likelihood maximization [14, 15]. The canonical way to use a mean- field 
approximation for the inverse problem is to calculate the correlations by linear response 
of the magnetizations mi to changes of the local fields hi to yield the connected 
correlations and solve for the couplings. In the same way, our approach is based on 
the Bethe-Peierls approximation [16] combined with the hnear response relation. 

We begin by deriving the connected correlation for a system of spins in a tree 
graph, where the Boltzmann distribution (1) can be written exactly in terms of the 
Bethe ansatz, that is the combination of the one-point marginals hi at the spins and the 
two-point marginals hij at the bonds, 

Pbp [s] = Wh^Sif-'^ n ' (4) 

where Zi is the number of non-zero bonds connected to spin i [16, 17, 18, 19, 20]. 
The marginals hi and hij, which are still unknown, can be parameterized by local 
magnetizations mj and correlation parameters Cij as 

= ^(1 + rUiSi) , 

hij{si, Sj) = ^[{1 + miSi){l + rUjSj) + CijSiSj] , (5) 
with constraints 

-1 < nii < +1 , 

— l + \mi + mj\— niimj < Cij < +1 — \mi — mj\ — niimj . (6) 

Instead of single spins coupled to effective fields in naive mean-field theory, the 
Bethe ansatz is based on bonds. It is thus a natural approximation to take in the 
context of the inverse problem, where it is bonds that are to be determined, not the 
statistics of spins. 

The self-consistent equations for the parameters and Cij describing the statistics 
of spin pairs in equation (5) can be found by minimizing the KuUback-Leibler divergence 
V between the Bethe ansatz (4) and the Boltzmann measure (1), X'(pBp[s],p[s]) = 
J2sPbp[s\ \n{pBp[s]/p[s]), giving 

hi+ JijTUj = (1 — Zj)arcth(mj) -|- 

Si + nijSiSj {1 + miSi){l + nijSj) + CijSiSj 

+ L — r — 4^ — ' 

7 "SjSj {I + miSi){l + mjSj) + CijSiSj 

Jio = L ^ In '—^ , (8) 

where di stands for the boundary set containing all spins that are connected to spin i 
by a bond. 
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The self-consistent equation (8) is quadratic in Cij, but only one of its solutions is 
compatible with the constraint (6), which is 

(1 - - m2)ti • + 2mimj 

Cij = '—1= mimj , (9) 

1 + ^D,j 

where Dij = 1 — 2mimjtij — (1 — mf — rn'j)tfj and tij = th(2Jjj) . 

The correlation parameters Cij give the connected correlation for pairs of spins 
that interact with each other in the tree. To find the connected correlation Cij between 
any pair of spins, we follow a route based on the linear response relation [4, 5, 20]. 
Differentiating the self-consistent equations (7) and (9) with respect to the fields hi 
yields a set of equations for the susceptibilities drui/dhj. The general linear response 
relation Cij = dnii/dhj then directly gives the connected correlations, for which we 
obtain 

C- 

{C-% = - J,, + 4- + ^^^^.)2_(l_'^2)(i_^2) A^^J) , (10) 

with 

J ^ 1 [ [(1 + m,)(l + m,) + Ci,][{l - m,)(l - m,) + Q,] \ 

4 \[(l+m,)(l-m,)-a,][(l-m,)(l+m,)-a,]j ' ^ ' 
where the coupling matrix Jij was extended to every pair of spins so that Jij = if z 
and j are not connected by a bond. 

Note that the first two terms in the right hand side of equation (10) cancel exactly 
due to (8), we however keep them to relate the expression to the Sessak-Monasson 
expansion later on. 

An equivalent version of equations (10)-(12) was derived by Welling and Teh in 
the context of the forward problem [5], and was used to estimate the data evidence 
in Bayesian inference [21]. These equations can be used to estimate the couplings and 
local fields in the inverse Ising problem. To do so, we first assume a fully connected 
model without self-interactions. It follows that Zi = N — 1, where is the number 
of spins. The Bethe ansatz, together with the equations (7), (8), (10) and (11), is 
exact only in tree graphs, using it as an approximation in loopy graphs is usually 
referred to as Bethe-Peierls approximation. In loopy graphs, the correlation between 
interacting spins estimated by the linear response relation (10) is different from the 
correlation parameter (9), which is, in general, of lower accuracy [G]. In equation (10), 
we identify the magnetizations mj and the connected correlations Cij with the values 



D 



and C^j estimated from data. Inserting the solution (9) for the correlation parameter 
Cij, equation (10) is solved for the couplings Jij. To calculate the local fields hi, the 
solutions Jij are fed into the self-consistent equation (7). Its solution in closed form 
completes the Bethe reconstruction. 

The Bethe-Peierls approximation is related to a number of previous approaches. 
Expanding equation (10) to first and second order in J^j recovers the naive mean-field 
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and TAP reconstructions, respectively. Calculating the couplings and local fields using 
the self-consistent equations (7) and (8) with = mf and Cij = estimated from 
data leads to the independent-pair approximation [s]. On the other hand, inserting 
the estimated values mf and into the places of mi and Cij as well as Cij in 
equation (10) (without using the solution (9) for Cij), one obtains the Sessak-Monasson 
reconstruction [8, 10]. And lastly, since the extrema of Bethe-Peierls free energy are the 
fixed points of belief propagation [22], susceptibility propagation is expected to have a 
fixed point at the solution of Bethe reconstruction. 

Interestingly, for a tree at zero field, the magnetizations nii are zero for finite system, 
and equations (10) and (11) take on a particularly simple form {C~^)ij = — sh( Jjj)ch( Jjj) 
for i 7^ j, and {C~^)ii = 1 + Z^jgOi sh^( Jjj) (see [23] for the discussion on symmetry 
breaking in Cayley trees in the thermodynamic limit). This result can be understood 
as a geometrical relation, where the connected correlation has a simple interpretation: 
For any two nodes i and j in a tree there is a unique path Pij connecting them, and 
the correlation function can be written as Cij = Yli^fzp-j th(Ji^) for i j, and Cu = 1, 
where denotes links in this path. Note that in the expression of Cij, if one replaces 
the multiplication by summation over the links, the correlation matrix then becomes 
the distance matrix in the weighted tree with weight Wij = th(Jjj). A similar expression 
for the inverse of the distance matrix is already well-known in graph theory [24]. Bethe 
reconstruction then turns out to be particularly simple, 

J,, = -iarcsh[2(C-i),,] , (z ^ j) . (13) 

It is clear from the above description that the computational expense of Bethe 
reconstruction is mostly due to the inversion of the correlation matrix, as is the case for 
MF, TAP and SM. These methods are low computational expense methods, which are 
convenient to apply to systems of large sizes. The expense of inversion of the correlation 
matrix is of the order of A^^ for methods such as Gaussian elimination. On the other 
hand, the computation in a single update of SusP is already of the order of A^^, which 
makes SusP slower due to multiple iterative steps. 

In the remainder of this note, we discuss the performance of Bethe reconstruction 
compared to MF, TAP, SM and SusP, which we refer to as the mean-field family. 
Although IP is faster than those methods, it usually leads to estimates of couplings 
inferior to all of these methods, so no comparison with IP is made. 

We will first compare the performance of these methods for a tree and for the 
Sherrington-Kirkpatrick (SK)-model [25] in the absence of sampling noise (effectively 
infinite number of patterns). Then the effect of sampling noise is considered for a 
particular case of sparse random graphs. The comparisons proceed as follows: We first 
construct a graph of A^ nodes, on which spins are located. Random values of couplings 
are assigned to each edge, and random local fields are assigned at each spin. 
Two distributions for the random variables, the normal distribution with zero mean 
and standard deviation a (A/'(0,cr^)) and the uniform distribution on the interval [a,b] 
{U{a,b)), will be used (see the corresponding figure captions). The magnetizations 
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mf and connected correlations are calculated either by enumerating all the spin 
configurations (infinite sampling) or by performing Monte-Carlo simulation of the model 
at inverse temperature /3 to generate a finite number of M samples (finite sampling). 
The observed magnetizations and the connected correlations then serve as inputs to 
estimate the model parameters (3hi and (3Jij (where we explicitly reintroduced the 
inverse temperature). Concentrating on the reconstruction of the couplings, we measure 
the deviation of each solution from the underling couplings by the relative deviation 



Figure 1 (a) shows the reconstructions of the Ising model on a Cayley tree of 
connectivity z = 3 with N = 22 spins. Here, the quality of Bethe reconstruction is 
limited only by the precision of numerical computations, confirming its exactness on a 
tree if the magnetizations and the connected correlations C^* are known exactly. 
SusP is also exact in this case, but is numerically less efficient and known to fail at low 
temperatures [12]. On the other hand MF, TAP, and SM are all approximate, as seen 
particularly when the couplings are strong. 

In figure 1 (b) we turn to loopy graphs, and consider a particularly extreme case, the 
fully connected SK-model of = 20 spins. At weak couplings (small /3), the difference 
between the reconstruction of all the methods and the underling couplings are small, 
and in practice the differences between the methods will be obscured by sampling noise 
(see below). On the other hand, at strong couplings (larger /3), Bethe reconstruction 
clearly outperforms all the other methods of the mean-field family. Interestingly, even 
inside the glassy regime, a finite overlap between the Bethe-reconstructed couplings 
and the underlying couplings persists. Note that SusP also follows the line of Bethe 
reconstruction closely at weak couplings, but at strong couplings it fails to converge to 
the Bethe solution. 

We now consider the Ising model with A^ = 50 spins placed on the vertices of 
a random graph with fixed connectivity z = 3 [26]. Such graphs typically contain 
loops of length In(A^) [27, 28], but retain a locally tree-like structure. Figure 2 (a) 
shows the results for a finite number of M = 5000 samples, where we omit SusP from 
the comparison due to its numerical instability. At weak couplings (corresponding 
to high temperatures), reconstruction by any method is limited by sampling noise. 
Remarkably, the sampling noise tends to smear out the small difference between SM 
and Bethe reconstructions for the entire weak coupling regime. Again, the deviation d 
of Bethe reconstruction increases only slowly through the strong coupling region, making 
it the best candidate for large couplings. Figure 2 (b) probes the performances of the 
different methods at /3 = 1.5 with varying number of samples. The results show that 
all the reconstruction methods from the mean-field family are sensitive to the sampling 
noise. For the present system, it requires more than 700 samples to obtain a good 
reconstruction and to see a clear difference between the different methods. 

In conclusion, we showed that the Bethe-Peierls approximation and linear response 
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can serve as the basis for a method of the mean-field family to solve the inverse 
Ising problem. This Bethe reconstruction is not only computationally efficient but 
also stable through a wide range of couplings. While MF, TAP can be considered as 
low-order expansions in the couplings of Bethe reconstruction, the Sessak-Monasson 
reconstruction can be recovered by treating the parametric correlations in a particular 
manner. Although our observations suggest that susceptibility propagation and Bethe 
reconstruction are closely related, further work is needed to probe their relationship, in 
particular in the regime of strong couplings. 

Since the Bethe-Peierls approximation works well for locally tree-like graphs, the 
extension of Bethe reconstruction to short-loop graphs, possibly made by linear response 
applied to the Kikuchi approximation [17, 19, 29], is an interesting direction for future 
research. Another direction is extending the method to non-binary degrees of freedom, 
for which the linear response relation has been already studied [11]. 
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Figure 1. Benchmarking in the infinite sampling limit. Performance of naive 
mean- field (MF), TAP (TAP) (the same as MF in these cases), Sessak-Monasson 
(SM), susceptibility propagation (SusP) and Bethe reconstruction (Bethe), measured 
by the relative deviation d (see main text) between the reconstructed couplings and 
the underlying couplings, (a) Cayley tree of connectivity z = 3 over N = 22 spins, 
J°- - W(-l, +1), h° = 0. (b) SK-model of ^ = 20 spins, J°- - Af{0,l/N), h° = 0. 




Figure 2. Benchmarking in the case of finite sampling. Performance of 
naive mean-field (MF), TAP (TAP), Sessak-Monasson (SM) and Bethe reconstruction 
(Bethe) on a random graph with fixed connectivity z = 3, iV = 50 spins, ~ 
U{—1,+1), ~ W(— 0.1, +0.1). (a) The deviation d versus the inverse temperature 
/?, with M = 5000 samples, (b) The deviation d versus the number of samples M at 
P = 1.5. 
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